% This code goes through various levels of BEA disaggregation, based on the
% actual inputs for price stickiness, sector size and I/O linkages
% by Raphael Schoenle, February 1, 2016 and revised April 2019

% note: one needs to run concording.m first

% only focus on case 6, our base case, across BEA aggregations

clearvars -except wd

% parameter sets: only base case for Figure 1!
%Parameters values
params.rho_a = 0.9;         % persistence agg tech shocks
params.rho_mu = 0.9;         % persistence monetary shock
params.rho_k = 0.9;         % persistence idiosync tech shock
params.beta = 0.9975;      % discount factor
params.theta = 6;          % elast of subst within sectors
params.delta = 0.5;

delta_set = [0.5];
eta_set = [2];
phi_set = [2];
sigma_set = [1];
phi_pi_set = [1.24];
phi_c_set = [0.33]./12;
phi_rho_i_set = [0];
phi_gc_i_set = [0];

n_parameters = length(phi_rho_i_set);
num_case = 18;
n_periods= 30;

% Position indices for shocks
pos_mushock = 1;

% produce aggregation
% set of disaggregation levels
kl = [ 341 56 7 ];

%% Sectors
for i_k=1:length(kl)
    kkk = kl(i_k);
    % set the number of sectors for the economy
    n_sectors = kkk;

    % load data
    clear FPA FPA_het Omega_C Omega Omega_new2 FPA_BEA_new Omega_C_new
	if n_sectors < 341
        filename = ['empirical_aggregation' num2str(n_sectors) '.mat']
        load(filename);

        FPA(:,:,6) = FPA_BEA_new;
        Omega_C(:,:,6) = Omega_C_new;
        Omega(:,:,6) = Omega_new2;
    else
        % 341 original most disaggregated data        
        
        Ck = csvread("2002CkDetailed.csv");
        Znum = csvread("2002RevShareDetailed.csv");
        n = length(Ck); % Number of Industries
        Z = Znum(1:n,1:n);

        Omega_C(:,:,6) = Ck/sum(Ck);

        % Data
        load FPA.mat
        FPA(:,:,6) = FPA;
        
        W = Z./(ones(n,1)*sum(Z)); %Normalize the share of inputs for each sector, omega
        Omega(:,:,6) = W;

    end

    % outer loop over the parameter values, index k

for k = 1:1 
    params.eta = eta_set(k);
    params.frisch = phi_set(k);
    params.sigma = sigma_set(k);
    params.phi_pi = phi_pi_set(k);
    params.phi_c = phi_c_set(k);
    params.phi_gc = phi_gc_i_set(k);
    params.rho_i = phi_rho_i_set(k);

    % inner loop over the cases, index i

    for i = 6
    i
    tic
    if i < num_case/3+1
        % n_sector cases
        %%Generate Gensys matrices
        [G0, G1, C, PSI, PI, pos_c, pos_pc, pos_mc, pos_yk, pos_p, pos_ps]  = ...
                Gensys_matrix(n_sectors, FPA(:,:,i),Omega_C(:,:,i), Omega(:,:,i), params);
    toc
        %%Call Gensys
        [G1, C, imp, fmat, fwt, ywt, gev, eu, loose] = gensys(G0, G1, C, PSI, PI);
    toc
        %%Generate IRFs
        %%%Monetary policy shocks
        % now include MC as output

            %check existence and uniqueness
            eu
        
            yt_imp_mu = Gensys_IRF(G1, n_sectors, n_periods, pos_mushock, imp);
            % save Gensys solution matrices
            filename = ['temp' num2str(i) '.mat'];
            save(filename,'yt_imp_mu');

            % output = [index for sectors, consumption IRF, inflation IRF, real MC IRF]
            output = [1:n_periods; yt_imp_mu(pos_c,2:end); ...
                yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1); ...
                yt_imp_mu(pos_mc,2:end) - (yt_imp_mu(pos_pc,2:end)) ]';
            autocorr_c = autocorr(output(:,end-1),1);
            autocorr_inf = autocorr(output(:,end),1);
            autocorr_mc = autocorr([mean(yt_imp_mu(pos_mc,2:end)) - (yt_imp_mu(pos_pc,2:end))]',1);

            % output matrix: initial C, Pi, MC, cumulative C, Pi, mean(MC)...
            % AR(1) of C, Pi, MC + IRFs (C, Pi, MC)
            output_matrix = [output(:,1) repmat([output(1,2:3)...
                mean(yt_imp_mu(pos_mc,2) - (yt_imp_mu(pos_pc,2)))...
            sum(output(:,2:3)) sum(yt_imp_mu(pos_mc,2:end) - (yt_imp_mu(pos_pc,2:end)))...    
            autocorr_c(2) autocorr_inf(2) autocorr_mc(2)],n_periods,1) output(:,2:end)];

            % save output
            filename = ['Gensys_mu_output_disaggregation_empirical_case_' num2str(i) 'params' num2str(k) '_' num2str(n_sectors) '.csv'];
            fid = fopen(filename, 'w');
            fprintf(fid,'Period\tConsumption_impact_irf_mu\tInflation_impact_irf_mu\tMC_mean_impact_irf_mu\tConsumption_cum_irf_mu\tInflation_cum_irf_mu\tMC_mean_cum_irf_mu\tConsumption_persistence_mu\tInflation_persistence_mu\tMC_persistence_mu\tConsumption_mu\tInflation_mu\tMC_mu \n');
            fclose(fid);
            dlmwrite(filename, (output_matrix), '-append', 'precision', '%.10f', 'delimiter', '\t');

            % Table 5 output 341 vs 56 sectors
            Table_5 = [sum(yt_imp_mu(pos_c,2:end)) ; ...
                sum(yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1)); ...
                sum(mean(yt_imp_mu(pos_mc,2:end) - (yt_imp_mu(pos_pc,2:end))))]
            filename = ['Table_5_' num2str(n_sectors) '.csv'];
            fid = fopen(filename, 'w');
            fclose(fid);
            dlmwrite(filename, (Table_5), '-append', 'precision', '%.10f', 'delimiter', '\t');
 
            % save output as matlab file for later plotting
            %filename = ['Gensys_mu_output_main_case_' num2str(i) 'params' num2str(k) '_341.mat'];
            %save(filename, 'output_matrix')
        
    end
    
    %%Save Case output in output matrix
    c_imp(:,i,1,kkk) = yt_imp_mu(pos_c,2:end)';
    inf_imp(:,i,1,kkk) = (yt_imp_mu(pos_pc,2:end) - yt_imp_mu(pos_pc,1:end-1))';
    mc_imp(:,i,1,kkk) = [mean(yt_imp_mu(pos_mc,2:end)) - (yt_imp_mu(pos_pc,2:end))]';
    
    toc
    end
end
%     
end
  

%     % make figures for N-sector economies combines
     i = 6
     figure()
     NameArray = {'LineStyle', 'Color','Marker'}
     ValueArray = {'-','b','none'; '--','b','.'; ':','b','.'; '-.', 'b','.' ;'-','b','s'; '--','b','+';'-.','b','d'}
    % consumption
     subplot(1,1,1)
     line = plot(1:n_periods,c_imp(:,i,1,kl(1))','linewidth',2)
     set(line ,NameArray, ValueArray(1,:))


     if length(kl)>1
        hold on

        for ll = 2:length(kl)
            line = plot(1:n_periods,c_imp(:,i,1,kl(ll))','linewidth',2)
            set(line, NameArray, ValueArray(ll,:))
        end
     end
     
     title('\textbf{Consumption}','Interpreter','latex')
     xlabel('\textbf{Period}','Interpreter','latex')
     ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
     dlegend=cell(1,length(kl));
    for hh=1:length(kl)
       dlegend(hh)={[ 'n=' num2str(kl(hh))]}
    end
    legend1=legend(dlegend,'Location','best')
    fig341 = gcf
    fig341.PaperUnits = 'normalized';
    fig341.PaperPosition = [0 0 1 .6];
    fig341.PaperPositionMode = 'manual';
    filename = ['IRFs_Nsectors_empirical_aggregation_C.pdf'];
    saveas(fig341,filename);
    
    % inflation
    figure()
    subplot(1,1,1)
    line = plot(1:n_periods,inf_imp(:,i,1,kl(1))','linewidth',2)
    set(line, NameArray, ValueArray(1,:))

     if length(kl)>1
        hold on

        for ll = 2:length(kl)
            line = plot(1:n_periods,inf_imp(:,i,1,kl(ll))','linewidth',2)
            set(line, NameArray, ValueArray(ll,:))
            end
     end
     
    title('\textbf{Inflation}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    dlegend=cell(1,length(kl));
    for hh=1:length(kl)
       dlegend(hh)={[ 'n=' num2str(kl(hh))]}
    end
    legend1=legend(dlegend,'Location','best')
    fig341 = gcf
    fig341.PaperUnits = 'normalized';
    fig341.PaperPosition = [0 0 1 .6];
    fig341.PaperPositionMode = 'manual';
    filename = ['IRFs_Nsectors_empirical_aggregation_infl.pdf'];
    saveas(fig341,filename);
    
    % marginal cost
    figure()
    subplot(1,1,1)
    line = plot(1:n_periods, mc_imp(:,i,1,kl(1))','linewidth',2)
    set(line, NameArray, ValueArray(1,:))
  
    if length(kl)>1
        hold on

        for ll = 2:length(kl)
            line = plot(1:n_periods, mc_imp(:,i,1,kl(ll))','linewidth',2)
            set(line, NameArray, ValueArray(ll,:))
        end
    end
    title('\textbf{Real MC}','Interpreter','latex')
    xlabel('\textbf{Period}','Interpreter','latex')
    ylabel('\textbf{Percentage Deviation from SS}','Interpreter','latex')
    dlegend=cell(1,length(kl));
    for hh=1:length(kl)
       dlegend(hh)={[ 'n=' num2str(kl(hh))]}
    end
    legend1=legend(dlegend,'Location','best')
    fig341 = gcf
    fig341.PaperUnits = 'normalized';
    fig341.PaperPosition = [0 0 1 .6];
    fig341.PaperPositionMode = 'manual';
    filename = ['IRFs_Nsectors_empirical_aggregation_MC.pdf'];
    saveas(fig341,filename);

save('alloutput_empirical_aggregation.mat')    


